rm(list=ls())
# Set the working directory to where folders named after the samples are located.
# The folder contains spliced.mtx, unspliced.mtx, barcodes and gene id files, and json files produced by scKB that documents the sequencing stats.
setwd("/scratch/AG_Ohler/CheWei/scKB")
# Load libraries
library(edgeR)
library(tidyverse)
library(Seurat)
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(GeneOverlap)
library(gprofiler2)
library(ggrepel)
library(ggplot2)
library(scater)
library(RColorBrewer)
library(future)
#for 300gb ram
options(future.globals.maxSize = 300000 * 1024^2)
rc.integrated <- readRDS("./rc.integrated_wt_shr_scr.rds")
feature_names <- read_tsv("./supp_data/features.tsv.gz", col_names = c("AGI", "Name", "Type")) %>%
select(-Type) %>%
distinct()
bscs <- read.csv("./Benfey_single_cell-Samples.csv", na.strings=c("","NA"), stringsAsFactors = F)
bscs <- bscs %>% select(c('sample','name','source','geno', 'genotype','transgene','treatment','age','timepoint','rep','target_cells','date','seq_run'))
bscs$date <- gsub('^([0-9]{4})([0-9]{2})([0-9]+)$', '\\1-\\2-\\3', bscs$date)
bscs$target_cells <- prettyNum(bscs$target_cells, big.mark = ',')
orig.idents <- data.frame(sample=rc.integrated$orig.ident)
bscs %>% filter(sample %in% orig.idents$sample)
orig.idents_info <- left_join(orig.idents, bscs, by="sample")
orig.idents_info
treatment <- orig.idents_info$treatment
names(treatment) <- colnames(rc.integrated)
head(treatment)
rc.integrated <- AddMetaData(
object = rc.integrated,
metadata = treatment,
col.name = 'treatment')
timepoint <- orig.idents_info$timepoint
names(timepoint) <- colnames(rc.integrated)
rc.integrated <- AddMetaData(
object = rc.integrated,
metadata = timepoint,
col.name = 'timepoint')
age <- orig.idents_info$age
names(age) <- colnames(rc.integrated)
rc.integrated <- AddMetaData(
object = rc.integrated,
metadata = age,
col.name = 'age')
rep <- orig.idents_info$rep
names(rep) <- colnames(rc.integrated)
rc.integrated <- AddMetaData(
object = rc.integrated,
metadata = rep,
col.name = 'rep')
date <- orig.idents_info$date
names(date) <- colnames(rc.integrated)
rc.integrated <- AddMetaData(
object = rc.integrated,
metadata = date,
col.name = 'date')
geno <- orig.idents_info$geno
names(geno) <- colnames(rc.integrated)
rc.integrated <- AddMetaData(
object = rc.integrated,
metadata = geno,
col.name = 'geno')
table(rc.integrated$treatment)
table(rc.integrated$geno)
rc.integrated$geno_trt <- rep("WT")
rc.integrated$geno_trt[rc.integrated$geno=="WT"] <- "WT"
rc.integrated$geno_trt[rc.integrated$geno=="shr-2"] <- "shr"
rc.integrated$geno_trt[rc.integrated$geno=="scr-4"] <- "scr"
rc.integrated$geno_trt <- factor(rc.integrated$geno_trt, levels=c("WT", "shr", "scr"))
table(rc.integrated$geno_trt)
table(rc.integrated$celltype.anno)
# Simple QC label
rc.integrated$celltype.anno <- gsub("Putative Quiescent Center", "Quiescent Center", rc.integrated$celltype.anno, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
order <- c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium", "Unknown")
palette <- c("#9400D3","#DCD0FF", "#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#FF9900","#E6194B", "#9A6324", "#FFE119","#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
options(repr.plot.width=6, repr.plot.height=6)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols=color)
options(repr.plot.width=12, repr.plot.height=6)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", split.by="geno_trt", cols=color)
sample_num_df <- data.frame(table(rc.integrated$geno_trt))
sample_num_df <- dplyr::rename(sample_num_df, "total_cells"="Freq", "orig.id"="Var1")
sample_num_df
rc.integrated$sample.cell.p <- paste(rc.integrated$geno_trt, rc.integrated$celltype.anno, sep = ",")
sample_cell_df <- data.frame(table(rc.integrated$sample.cell.p))
sample_cell_df <- separate(sample_cell_df, col="Var1", sep=",", into=c("orig.id", "celltype")) %>%
dplyr::rename("n_celltype"="Freq") %>%
left_join(sample_num_df) %>%
mutate(proportion=n_celltype/total_cells)
sample_cell_df
wt_cells <- filter(sample_cell_df, orig.id == "WT") %>%
group_by(celltype) %>%
summarise(wt_prop=mean(proportion))
wt_cells
sample_cell_df <- left_join(sample_cell_df, wt_cells) %>%
mutate(wt_norm_prop=proportion/wt_prop)
sample_cell_df
sample_cell_wide <- sample_cell_df %>%
ungroup() %>%
select(orig.id, celltype, wt_norm_prop) %>%
pivot_wider(names_from = celltype, values_from = wt_norm_prop)
sample_cell_wide
sample_cell_wide[is.na(sample_cell_wide)==T] <- 0.0000001
sample_cell_wide
cell_m <- as.matrix(log2(sample_cell_wide[, -1]))
rownames(cell_m) <- sample_cell_wide$orig.id
cell_m
# add to hm cell_fun = function(j, i, x, y, width, height, fill) {grid.text(sprintf("%.0f", int_matrix[i, j]), x, y, gp = gpar(fontsize = 10))})
Cell_hm <- Heatmap(cell_m, name = "log2(mutant/WT cell proportion)",
heatmap_legend_param = list(title_position="topcenter",
color_bar = "continuous",
legend_direction = "horizontal"),
col = colorRamp2(c(-3, 0, 3), c("blue", "beige", "#e31a1c")),
cluster_rows = F,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE, show_row_names = TRUE, show_row_dend = TRUE, show_column_dend = TRUE,
clustering_distance_rows = "pearson", clustering_distance_columns = "pearson",
row_names_gp = gpar(fontsize = 12))
draw(Cell_hm, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top")
sample_num_df <- data.frame(table(rc.integrated$orig.ident))
sample_num_df <- dplyr::rename(sample_num_df, "total_cells"="Freq", "orig.id"="Var1")
sample_num_df
rc.integrated$sample.cell.p <- paste(rc.integrated$orig.ident, rc.integrated$celltype.anno, sep = ",")
sample_cell_df <- data.frame(table(rc.integrated$sample.cell.p))
sample_cell_df <- separate(sample_cell_df, col="Var1", sep=",", into=c("orig.id", "celltype")) %>%
dplyr::rename("n_celltype"="Freq") %>%
left_join(sample_num_df) %>%
mutate(proportion=n_celltype/total_cells)
sample_cell_df
wt_cells <- filter(sample_cell_df, orig.id %in% c("sc_12", "sc_20", "sc_21", "sc_30", "sc_31", "sc_51")) %>%
group_by(celltype) %>%
summarise(wt_prop=mean(proportion))
wt_cells
sample_cell_df <- left_join(sample_cell_df, wt_cells) %>%
mutate(wt_norm_prop=proportion/wt_prop)
sample_cell_df
sample_cell_wide <- sample_cell_df %>%
ungroup() %>%
select(orig.id, celltype, wt_norm_prop) %>%
pivot_wider(names_from = celltype, values_from = wt_norm_prop)
sample_cell_wide
sample_cell_wide[is.na(sample_cell_wide)==T] <- 0.0000001
sample_cell_wide
cell_m <- as.matrix(log2(sample_cell_wide[, -1]))
rownames(cell_m) <- sample_cell_wide$orig.id
cell_m
options(repr.plot.width=12, repr.plot.height=10)
# add to hm cell_fun = function(j, i, x, y, width, height, fill) {grid.text(sprintf("%.0f", int_matrix[i, j]), x, y, gp = gpar(fontsize = 10))})
Cell_hm <- Heatmap(cell_m, name = "log2(sample/WT cell proportion)",
heatmap_legend_param = list(title_position="topcenter",
color_bar = "continuous",
legend_direction = "horizontal"),
col = colorRamp2(c(-2, 0, 2), c("blue", "beige", "#e31a1c")),
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE, show_row_names = TRUE, show_row_dend = TRUE, show_column_dend = TRUE,
clustering_distance_rows = "pearson", clustering_distance_columns = "pearson",
row_names_gp = gpar(fontsize = 12))
draw(Cell_hm, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top")
table(rc.integrated$batch)
table(rc.integrated$batch, rc.integrated$orig.ident)
rc.integrated$batch<- gsub(" ", "_", rc.integrated$batch, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
table(rc.integrated$batch, rc.integrated$orig.ident)
# subset only samples you want to compare
integrated.de <- subset(rc.integrated, subset = geno_trt %in% c("WT", "shr"))
DefaultAssay(integrated.de) <- "RNA"
merged <- as.SingleCellExperiment(integrated.de)
merged
abundances <- table(merged$celltype.anno, merged$orig.ident)
abundances <- unclass(abundances)
abundances
# Attaching some column metadata.
extra.info <- colData(merged)[match(colnames(abundances), merged$orig.ident),]
y.ab <- DGEList(abundances, samples=extra.info)
y.ab
keep <- filterByExpr(y.ab, group=y.ab$samples$orig.ident)
y.ab <- y.ab[keep,]
summary(keep)
(design <- model.matrix(~batch + factor(geno_trt), y.ab$samples))
y.ab <- estimateDisp(y.ab, design, trend="none")
summary(y.ab$common.dispersion)
plotBCV(y.ab, cex=1)
fit.ab <- glmQLFit(y.ab, design, robust=TRUE, abundance.trend=FALSE)
summary(fit.ab$var.prior)
summary(fit.ab$df.prior)
plotQLDisp(fit.ab, cex=1)
res <- glmQLFTest(fit.ab, coef=ncol(design))
summary(decideTests(res))
topTags(res, n = nrow(res$table))
shr_DA_results <- topTags(res, n = nrow(res$table))$table
shr_DA_results
# subset only samples you want to compare
integrated.de <- subset(rc.integrated, subset = geno_trt %in% c("WT", "scr"))
DefaultAssay(integrated.de) <- "RNA"
merged <- as.SingleCellExperiment(integrated.de)
merged
abundances <- table(merged$celltype.anno, merged$orig.ident)
abundances <- unclass(abundances)
abundances
# Attaching some column metadata.
extra.info <- colData(merged)[match(colnames(abundances), merged$orig.ident),]
y.ab <- DGEList(abundances, samples=extra.info)
y.ab
keep <- filterByExpr(y.ab, group=y.ab$samples$orig.ident)
y.ab <- y.ab[keep,]
summary(keep)
(design <- model.matrix(~batch + factor(geno_trt), y.ab$samples))
y.ab <- estimateDisp(y.ab, design, trend="none")
summary(y.ab$common.dispersion)
plotBCV(y.ab, cex=1)
fit.ab <- glmQLFit(y.ab, design, robust=TRUE, abundance.trend=FALSE)
summary(fit.ab$var.prior)
summary(fit.ab$df.prior)
plotQLDisp(fit.ab, cex=1)
res <- glmQLFTest(fit.ab, coef=ncol(design))
summary(decideTests(res))
topTags(res, n = nrow(res$table))
(scr_DA_results <- topTags(res, n = nrow(res$table))$table)
y.ab2 <- calcNormFactors(y.ab)
y.ab2$samples$norm.factors
y.ab2 <- estimateDisp(y.ab2, design, trend="none")
fit.ab2 <- glmQLFit(y.ab2, design, robust=TRUE, abundance.trend=FALSE)
res2 <- glmQLFTest(fit.ab2, coef=ncol(design))
summary(decideTests(res2))
topTags(res, n = nrow(res$table))$table
(scr_DA_results_adj <- topTags(res, n = nrow(res$table))$table)
(shr_DA_results <- rownames_to_column(shr_DA_results, var = "celltype"))
shr_DA_results$geno <- rep("shr")
(scr_DA_results <- rownames_to_column(scr_DA_results, var = "celltype"))
scr_DA_results$geno <- rep("scr")
(DA_results <- bind_rows(shr_DA_results, scr_DA_results))
(DA_results <- mutate(DA_results, star = case_when(FDR <=0.001 ~ "***",
FDR <=0.01 ~ "**",
FDR <=0.05 ~ "*",
TRUE ~ " ")))
DA_results
write.csv(DA_results, "./supp_data/shr_scr_DA_celltype.csv", row.names=F)
FC_wide <- select(DA_results, geno, , celltype, logFC) %>%
pivot_wider(names_from = geno, values_from = logFC) %>% mutate(WT=rep(0))
FC_wide_m <- as.matrix(FC_wide[, 2:ncol(FC_wide)])
rownames(FC_wide_m) <- FC_wide$celltype
FC_wide_m
star.wide <- select(DA_results, geno, , celltype, star) %>%
pivot_wider(names_from = geno, values_from = star) %>%
mutate(WT=rep(" "))
star.wide[is.na(star.wide)] <- " "
star.wide_m <- as.matrix(star.wide[, 2:ncol(star.wide)])
rownames(star.wide_m) <- star.wide$celltype
FC_wide_m
star.wide_m
options(repr.plot.width=6, repr.plot.height=12)
vert_Cell_hm <- Heatmap(FC_wide_m,
name = "Mutant/WT Cell type",
col = colorRamp2(c(-4, 0, 4), c("blue", "beige", "#e31a1c")),
column_title = " ",
cluster_rows = T,
cluster_columns = F,
clustering_distance_columns = "pearson",
clustering_distance_rows = "pearson",
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
show_column_dend = TRUE,
heatmap_legend_param = list(title_position="topcenter",
color_bar = "continuous",
legend_direction = "horizontal"),
cell_fun = function(j, i, x, y, width, height, fill)
{grid.text(sprintf("%s", star.wide_m[i, j]), x, y, gp = gpar(fontsize = 10))})
draw(vert_Cell_hm, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top")
options(repr.plot.width=3.6, repr.plot.height=9)
FC_wide_m_2 <- FC_wide_m[c("Endodermis",
"Cortex",
"Pericycle",
"Phloem",
"Xylem",
"Procambium",
"Quiescent Center",
"Trichoblast",
"Atrichoblast",
"Stem Cell Niche",
"Columella",
"Lateral Root Cap") ,
c("shr", "scr")]
star.wide_m_2 <- star.wide_m[c("Endodermis",
"Cortex",
"Pericycle",
"Phloem",
"Xylem",
"Procambium",
"Quiescent Center",
"Trichoblast",
"Atrichoblast",
"Stem Cell Niche",
"Columella",
"Lateral Root Cap") ,
c("shr", "scr")]
col_lab_it <- c(
expression(italic("shr ")),
expression(italic("scr ")))
vert_Cell_hm2 <- Heatmap(FC_wide_m_2,
name = "Cell Type logFC \n (mutant/WT)",
col = colorRamp2(c(-3, 0, 3), c("blue", "beige", "#e31a1c")),
column_title = " ",
cluster_rows = F,
cluster_columns = F,
clustering_distance_columns = "pearson",
clustering_distance_rows = "pearson",
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
show_column_dend = TRUE,
column_labels=col_lab_it,
column_names_gp = gpar(fontsize = 16),
row_names_gp = gpar(fontsize = 16),
heatmap_legend_param = list(title_position="topcenter",
color_bar = "continuous",
legend_direction = "horizontal",
labels_gp = gpar(fontsize = 16),
title_gp = gpar(fontsize = 16, font = 2),
legend_height = unit(4, "cm")),
cell_fun = function(j, i, x, y, width, height, fill)
{grid.text(sprintf("%s", star.wide_m_2[i, j]), x, y, gp = gpar(fontsize = 14))})
draw(vert_Cell_hm2, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top")
hm_final <- grid.grabExpr(draw(vert_Cell_hm2, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top"))
FC_wide_m_2 <- t(FC_wide_m_2)
rownames(FC_wide_m_2) <- c("shr-2", "scr-4")
star.wide_m_2 <- t(star.wide_m_2)
rownames(star.wide_m_2) <- c("shr-2", "scr-4")
FC_wide_m_2
star.wide_m_2
options(repr.plot.width=8, repr.plot.height=4)
hor_Cell_hm2 <- Heatmap(FC_wide_m_2,
name = "Cell Type logFC \n (mutant/WT)",
col = colorRamp2(c(-3, 0, 3), c("blue", "beige", "#e31a1c")),
column_title = " ",
cluster_rows = F,
cluster_columns = F,
clustering_distance_columns = "pearson",
clustering_distance_rows = "pearson",
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
row_names_side = "left",
show_column_dend = TRUE,
column_names_rot = 45,
heatmap_legend_param = list(title_position="topcenter",
color_bar = "continuous",
legend_direction = "horizontal"),
row_names_gp = gpar(fontsize = 14, font = 3),
cell_fun = function(j, i, x, y, width, height, fill)
{grid.text(sprintf("%s", star.wide_m_2[i, j]), x, y, gp = gpar(fontsize = 14))})
draw(hor_Cell_hm2, padding = unit(c(10, 20, 10, 10), "mm"), heatmap_legend_side = "top")
pdf("./supp_data/shr_scr_horizontal_DA.pdf", width = 8, height = 4)
draw(hor_Cell_hm2, padding = unit(c(10, 20, 10, 10), "mm"), heatmap_legend_side = "top")
dev.off()
anno_scores <- data.frame(geno=rc.integrated$geno_trt,
prediction.score.Endodermis=rc.integrated$prediction.score.Endodermis,
prediction.score.Cortex=rc.integrated$prediction.score.Cortex,
annotation=rc.integrated$celltype.anno,
timezone=rc.integrated$time.anno,
consensus.time=rc.integrated$consensus.time.group)
anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F)
pred_plot <- anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F) %>%
ggplot(aes(x=prediction.score.Endodermis, y=prediction.score.Cortex, color=annotation)) + geom_point(alpha=0.1)
library(ggridges)
options(repr.plot.width=8, repr.plot.height=6)
pred_plot <- anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F) %>%
mutate(Endodermis_Cortex_diff=prediction.score.Cortex-prediction.score.Endodermis) %>%
ggplot(aes(y=factor(geno, levels=c("scr", "shr", "WT")), x=Endodermis_Cortex_diff, fill=geno)) + geom_density_ridges(rel_min_height = 0.00)
mylabels <- c(
expression(italic("scr-4")),
expression(italic("shr-2")), "WT")
pred_plot <- pred_plot +
theme_cowplot(font_size = 20) +
xlab("Cortex Prediction score - Endodermis Prediction score") +
ylab("") +
scale_fill_viridis_d(labels = rev(mylabels)) +
scale_x_continuous(limits=c(-1.2, 1.2)) +
scale_y_discrete(labels = mylabels)
pred_plot
options(repr.plot.width=8, repr.plot.height=10)
pred_plot <- pred_plot + facet_wrap(~ factor(timezone, levels=c("Maturation","Elongation","Meristem")), ncol = 1) + theme(legend.title = element_blank())
pred_plot
# make triangles
options(repr.plot.width=8, repr.plot.height=2)
d=data.frame(x=c(1,1,16, 16,16,1), y=c(1,2,1.5, 1.7,2.7,2.2), t=c('a', 'a', 'a', 'b', 'b', 'b'))
d
triangle_label <- ggplot() +
geom_polygon(data=d, mapping=aes(x=x, y=y, group=t, fill=t)) + annotate("text", x = 3.5, y = 1.5, label = "Endodermis", size = 8, color="white") +
annotate("text", x = 14.5, y = 2.2, label = "Cortex", size = 8, color="white") +
scale_fill_manual(values=(c("#0000FF", "#82B6FF"))) + theme_nothing() + scale_x_continuous(limits=c(0.1, 17))
triangle_label
options(repr.plot.width=9, repr.plot.height=11)
(pred_tri <- plot_grid(pred_plot, triangle_label, ncol=1, rel_heights = c(1, .2)))
ggsave("./supp_data/WT_shr_scr_prediction_ridgeplot.pdf", width=9, height=10)
ggsave("./supp_data/WT_shr_scr_prediction_ridgeplot_taller.pdf", width=9, height=11)
ggsave("./supp_data/WT_shr_scr_prediction_ridgeplot_tallest.pdf", width=9, height=12)
options(repr.plot.width=15, repr.plot.height=6)
time_plt <- DimPlot(rc.integrated,
group.by = "time.anno",
order = c("Meristem","Elongation","Maturation"),
cols = c("#1A237E", "#42B3D5", "#DCEDC8"), split.by="geno_trt")
time_plt
options(repr.plot.width=12, repr.plot.height=12)
cell_plt <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", split.by="geno_trt", cols=color)
(cell_time_plt <- plot_grid(cell_plt, time_plt, nrow=2))
options(repr.plot.width=20, repr.plot.height=12)
plot_grid(cell_time_plt, pred_tri, ncol=2, rel_widths = c(2,1.2))
# swtich direction
options(repr.plot.width=8, repr.plot.height=6)
pred_plot <- anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F) %>%
mutate(Cortex__Endo_diff=prediction.score.Endodermis-prediction.score.Cortex) %>%
ggplot(aes(y=factor(geno, levels=c("shr", "scr", "WT")), x=Cortex__Endo_diff, fill=geno)) + geom_density_ridges(rel_min_height = 0.00)
mylabels <- c(
expression(italic("shr")),
expression(italic("scr")), "WT")
pred_plot <- pred_plot +
theme_cowplot() +
xlab("Endodermis Prediction score - Cortex Prediction score") +
ylab("") +
scale_fill_viridis_d(labels = rev(mylabels)) +
scale_x_continuous(limits=c(-1.2, 1.2)) +
scale_y_discrete(labels = mylabels) + theme(legend.position = "none")
pred_plot
options(repr.plot.width=8, repr.plot.height=10)
pred_plot <- pred_plot + facet_wrap(~ factor(timezone, levels=c("Maturation","Elongation","Meristem")), ncol = 1) + theme(legend.title = element_blank())
pred_plot
# make triangles
options(repr.plot.width=8, repr.plot.height=2)
d=data.frame(x=c(1,1,16, 16,16,1), y=c(1,2,1.5, 1.7,2.7,2.2), t=c('a', 'a', 'a', 'b', 'b', 'b'))
d
triangle_label <- ggplot() +
geom_polygon(data=d, mapping=aes(x=x, y=y, group=t, fill=t)) + annotate("text", x = 2.5, y = 1.5, label = "Cortex", size = 7, color="white") +
annotate("text", x = 13.8, y = 2.2, label = "Endodermis", size = 7, color="white") +
scale_fill_manual(values=(c("blue", "red"))) + theme_nothing() + scale_x_continuous(limits=c(0.1, 17))
triangle_label
options(repr.plot.width=8, repr.plot.height=10)
(pred_tri <- plot_grid(pred_plot, triangle_label, ncol=1, rel_heights = c(1, .15), align="v"))
rc.integrated$Endodermis_Cortex_diff <- rc.integrated$prediction.score.Cortex - rc.integrated$prediction.score.Endodermis
options(repr.plot.width=15, repr.plot.height=6)
FeaturePlot(rc.integrated, features = "Endodermis_Cortex_diff",
cols = c("blue", "grey", "red"), label=F, repel=F, pt.size = 0.5, split.by="geno_trt")
options(repr.plot.width=15, repr.plot.height=18)
# Visualize co-expression of two features simultaneously
FeaturePlot(rc.integrated, features = c("prediction.score.Cortex", "prediction.score.Endodermis"), blend = TRUE, split.by="geno_trt", cols = c("grey","blue", "red"))
WT_ground <- readRDS("./supp_data/Ground_Tissue_Atlas.rds")
scr_ground <- readRDS("./supp_data/Ground_Tissue_scr.rds")
shr_ground <- readRDS("./supp_data/Ground_Tissue_shr.rds")
options(repr.plot.width=15, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
FeaturePlot(scr_ground, features = c("prediction.score.Cortex", "prediction.score.Endodermis"), blend = TRUE, cols = c("grey","blue", "red"), min.cutoff = 0, max.cutoff = 1)
summary(scr_ground$prediction.score.Cortex)
options(repr.plot.width=15, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
FeaturePlot(shr_ground, features = c("prediction.score.Cortex", "prediction.score.Endodermis"), blend = TRUE, cols = c("grey","blue", "red"), min.cutoff = 0, max.cutoff = 1)
options(repr.plot.width=6, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
(scr_e <- FeaturePlot(scr_ground, features = c("prediction.score.Endodermis"), cols = c("grey", "red"), min.cutoff = 0, max.cutoff = 1) + ggtitle("scr \n Prediction Score \n Endodermis"))
options(repr.plot.width=6, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
(scr_c <- FeaturePlot(scr_ground, features = c("prediction.score.Cortex"), cols = c("grey", "blue"), min.cutoff = 0, max.cutoff = 1) + ggtitle("scr \n Prediction Score \n Cortex"))
scr_time_plt <- DimPlot(scr_ground,
group.by = "time.anno",
order = c("Meristem","Elongation","Maturation"),
cols = c("#1A237E", "#42B3D5", "#DCEDC8")) + ggtitle("scr \n Developmental Stage") + theme(plot.title = element_text(hjust = 0.5))
scr_time_plt
options(repr.plot.width=6, repr.plot.height=8)
scr_ct_plt <-DimPlot(scr_ground, reduction = "umap", group.by = "consensus.time.group", cols=brewer.pal(10,"Spectral")) + ggtitle("scr \n Consensus Time") + theme(plot.title = element_text(hjust = 0.5))
scr_ct_plt
options(repr.plot.width=6, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
(shr_e <- FeaturePlot(shr_ground, features = c("prediction.score.Endodermis"), cols = c("grey", "red"), min.cutoff = 0, max.cutoff = 1) + ggtitle("shr \n Prediction Score \n Endodermis"))
options(repr.plot.width=6, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
(shr_c <- FeaturePlot(shr_ground, features = c("prediction.score.Cortex"), cols = c("grey", "blue"), min.cutoff = 0, max.cutoff = 1) + ggtitle("shr \n Prediction Score \n Cortex"))
shr_time_plt <- DimPlot(shr_ground,
group.by = "time.anno",
order = c("Meristem","Elongation","Maturation"),
cols = c("#1A237E", "#42B3D5", "#DCEDC8")) + ggtitle("shr \n Developmental Stage") + theme(plot.title = element_text(hjust = 0.5))
shr_time_plt
options(repr.plot.width=6, repr.plot.height=8)
shr_ct_plt <-DimPlot(shr_ground, reduction = "umap", group.by = "consensus.time.group", cols=brewer.pal(10,"Spectral")) + ggtitle("shr \n Consensus Time") + theme(plot.title = element_text(hjust = 0.5))
shr_ct_plt
options(repr.plot.width=18, repr.plot.height=16)
(pred_dim_plots <- plot_grid(scr_time_plt, scr_c, scr_e, shr_time_plt, shr_c, shr_e, ncol=3, align="hv"))
options(repr.plot.width=24, repr.plot.height=12)
plot_grid(pred_dim_plots, pred_tri, ncol=2, rel_widths = c(2.5,1.2))
ggsave("./supp_data/scr_shr_pred_dim.pdf", width=30, height=12)
options(repr.plot.width=8, repr.plot.height=6)
pred_plot <- anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F) %>%
mutate(Cortex__Endo_diff=prediction.score.Endodermis-prediction.score.Cortex) %>%
ggplot(aes(y=factor(geno, levels=c("shr", "scr", "WT")), x=Cortex__Endo_diff, fill=geno)) + geom_density_ridges(rel_min_height = 0.00)
mylabels <- c(
expression(italic("shr")),
expression(italic("scr")), "WT")
pred_plot <- pred_plot +
theme_cowplot() +
xlab("Endodermis Prediction score - Cortex Prediction score") +
ylab("") +
scale_fill_viridis_d(labels = rev(mylabels)) +
scale_x_continuous(limits=c(-1.2, 1.2)) +
scale_y_discrete(labels = mylabels) + theme(legend.position = "none")
pred_plot
options(repr.plot.width=8, repr.plot.height=15)
con_time_order <- c("T0", "T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9")
con_time_order <- rev(con_time_order)
pred_plot <- pred_plot + facet_wrap(~ factor(consensus.time, levels=con_time_order), ncol = 1) + theme(legend.title = element_blank())
pred_plot
options(repr.plot.width=8, repr.plot.height=20)
(pred_tri <- plot_grid(pred_plot, triangle_label, ncol=1, rel_heights = c(1, .1)))
options(repr.plot.width=24, repr.plot.height=16)
(pred_dim_plots <- plot_grid(scr_time_plt, scr_ct_plt, scr_c, scr_e, shr_time_plt, shr_ct_plt, shr_c, shr_e, ncol=4, align="hv"))
options(repr.plot.width=30, repr.plot.height=15)
plot_grid(pred_dim_plots, pred_tri, ncol=2, rel_widths = c(3.5,1.2))
ggsave("./supp_data/scr_shr_pred_dim_with_consensus_time.pdf", width=30, height=15)